﻿//******************************************************************
//  caja.h
//  Objeto 'caja de simulacion'  para MDDipolar.cpp
//	Grupo de Fisica de Fluidos y Biocoloides - Univers. de Granada
//  ©2011, Dr. Juan Carlos Fernandez Toledano
//******************************************************************
class caja
{
private:

public:
	double lbox[3];//Lado de la caja
	particula *part; //Array de particulas
	int npart;//Numero de particulas
	double phi;//Concentracion en volumen
	double rho;//Concentracion en masa
	double temp;//Temperatura de caja
	double temp0;//Temperatura de referencia
	double vol;//Volumen
	double rcut;//Distancia de corte de interaccion
	double dr[3];//Distancia a pares entre particulas
	double drnb[3];//Distancia a pares entre particulas sin condiciones periodicas
	double presion;

	///Parametros lara la link-cell (vecinos)
	int ncell;//Numero de celdas 
	int ncell2;//ncell^2
	int ncell3;//ncell^3
	double lgrid;//Lado de celda
	/*head[i] contiene el primer indice de particula contenido en celda i
	head[i]=-1 si no hay particulas en celda i*/
	int *head;
	//Contiene la lista de particulas enlazadas
	int *link;
	/*Ejemplo:
	celda i
	p1=head[i]
		- p1 es el indice de la primera particula en celda i. 
		- si no hay particulas en celda i, p1=-1
	p2 = link[p1]
		- p2 es el indice de la segunda particula en celda i
		- p2=-1 si no hay mas particulas en la celda
	etc...
	Por lo tanto, en la celda i habrá
	[p1,p2,..., pk] particulas	
	*/
	//Energias
	double epoten;//Energia potencia por particula en la caja
	double virial;
	double ekin;//Energia cinetica por particula en la caja
	//Control de la velocidad (termostato)
	double alpha; //Coeficiente de transformacion de velocidades (sqrt(Tdeseada/Treal))
	int nacc;///Paso de la simulacion
	double vel[3];///Velocidad media de las particulas de la caja
	//Elementos del 'microscopic stress tensor' asociado a la componente cinetica
	double svxy, svxz, svyz, svzz;
	//Elementos del 'microscopic stress tensor' asociado al potencial de interaccion
	double spxy, spxz, spyz, spzz;
	//Campo magnetico
	double u0, u30;
	//Parametros de la simulacion
	parametros param;

	caja()
	{

	}

	//caja(int np, double xT, double xphi, double xu, double xrcut, double Dt)
	caja(parametros par)
	{
		///Empleamos paso de tiempo Dt para generar las posiciones posOld
		///de las particulas empleando las velocidades
		param=par;
		if(param.rest==0)
		{
			//Iniciamos la simulacion desde el principio
			u0=param.u0;
			u30=3.0*u0;
			npart=param.npart;
			temp0 =1.0;
			phi=param.phi;
			rcut = param.rcut;
			rho = phi*6.0/pi;
			double bx=pow((4./3.)*pi*npart/phi,1./3.);

			for(int i=0;i<3;i++)
				lbox[i]= bx;
			part = new particula[npart];
			link = new int[npart];
			SetMeshPart2();
			initForces();
			initVelocities(param.Dt);
			vol=lbox[0]*lbox[1]*lbox[2];

			ncell=int(lbox[0]/rcut);
			ncell2=ncell*ncell;
			ncell3=ncell*ncell*ncell;
			lgrid= lbox[0]/double(ncell);
			cout<<"Numero de celdas: "<<ncell<<endl;
			cout<<"Lado de celda: "<<lgrid<<endl;
			head = new int[ncell*ncell*ncell];
			alpha=1.0;
		}
		else
		{
			//Reiniciamos una simulacion desde archivo existente
			loadConfig();
		}
	}

	///Reinicia la simulacion desde un archivo
	caja(char filename[50], double Dt)
	{

		///Empleamos paso de tiempo Dt para generar las posiciones posOld
		///de las particulas empleando las velocidades
		ifstream MiArchivo (filename,ios::in);
		if (MiArchivo.is_open())
		{
			///Cargamos las caracteristicas de la caja
			MiArchivo>>nacc>>npart>>phi>>temp>>rcut>>lbox[0]>>lbox[1]>>lbox[2];
			rho = phi*6.0/pi;
			part = new particula[npart];
			link = new int[npart];
			SetMeshPart();
			initForces();
			initVelocities(Dt);
			vol=lbox[0]*lbox[1]*lbox[2];

			ncell=int(lbox[0]/rcut);
			ncell2=ncell*ncell;
			ncell3=ncell*ncell*ncell;
			lgrid= lbox[0]/double(ncell);
			cout<<"Numero de celdas: "<<ncell<<endl;
			cout<<"Lado de celda: "<<lgrid<<endl;
			head = new int[ncell*ncell*ncell];
			alpha=1.0;
			///Cargamos las posiciones de las particulas
			for(int i=0;i<npart;i++)
				for(int j=0;j<3;j++)
				{
					MiArchivo>>part[i].pos[j];
				}
				///Cargamos las velocidades
				for(int i=0;i<npart;i++)
					for(int j=0;j<3;j++)
					{
						MiArchivo>>part[i].vel[j];
						///Construimos las posiciones anteriores de las particulas
						part[i].posOld[j]=part[i].pos[j]-part[i].vel[j]*Dt;
					}
					MiArchivo.close();
		}

	}

	~caja()
	{
		delete [] part;
	}

	///Imprime la configuracion
	void printConfig(char namefile[30]);
	///Coloca las particulas de forma aleatoria en una caja
	void SetRandomPart(void);
	///Coloca las particulas en red
	void SetMeshPart(void);
	///Inicializa las fuerzas
	void initForces(void);
	///Imprime la informacion de la particula considerada
	void followPart(int ind);
	///Calcula la linkCell
	void linkCell(void);
	///Calcula la lista de vecinos
	void listaVecinos(void);
	///Inicializa fuerzasK suma ewald yukawa (No probado)
	void initForcesKR(void);
	///Inicializa las velocidades y calcula posOld
	void initVelocities(double Dt);
	///Aplica condiciones contorno periodicas a particula ind
	void boundary(int ind);
	///Calcula la distancia cuadratica entre particulas
	double distancia2(int ind1,int ind2);
	///Calcula la distancia entre particulas
	double distancia(int ind1,int ind2);
	///Imprime el archivo XYZ de la configuracion de particulas
	void printXYZfile(char namefile[30]);
	//Se imprime el tensor de esfuerzos microscopicos
	void printShearStressTensor(char namefile[30]);
	//Inicializacion de las velocidades
	void iniVelocidades2(void);
	void SetMeshPart2(void);
	void loadConfig(void);
};

void caja::loadConfig(void)
{
	ifstream MiArchivo (param.iniFile,ios::in);
	char line[150];
	char dummy[20];
	if(MiArchivo.fail())
	{
		cout<<"No se pudo encontrar el archivo de reinicio "<<param.iniFile<<endl;
		exit(1);
	}
	if (MiArchivo.is_open())
	{
		MiArchivo.getline(line,150);
		MiArchivo.getline(line,150);
		MiArchivo.getline(line,150);
		MiArchivo.getline(line,150);

//		MiArchivo.getline(line,150);


		MiArchivo>>param.npart>>param.phi>>param.rcut>>param.Dt>>param.rest>>param.ncorr
			>>param.trelaj1>>lbox[0]>>lbox[1]>>lbox[2]>>dummy;
		MiArchivo.getline(line,150);
		MiArchivo.getline(line,150);
		MiArchivo.getline(line,150);
		MiArchivo.getline(line,150);
		MiArchivo.getline(line,150);
		MiArchivo.getline(line,150);

		//Construimos las cosas
		u0=param.u0;
		u30=3.0*u0;
		npart=param.npart;
		temp0 =1.0;
		phi=param.phi;
		rcut = param.rcut;
		rho = phi*6.0/pi;
		part = new particula[npart];
		link = new int[npart];
		initForces();
		//double bx=pow((4./3.)*pi*npart/phi,1./3.);

		//for(int i=0;i<3;i++)
		//	lbox[i]= bx;
		vol=lbox[0]*lbox[1]*lbox[2];
		ncell=int(lbox[0]/rcut);
		ncell2=ncell*ncell;
		ncell3=ncell*ncell*ncell;
		lgrid= lbox[0]/double(ncell);
		cout<<"Numero de celdas: "<<ncell<<endl;
		cout<<"Lado de celda: "<<lgrid<<endl;
		head = new int[ncell*ncell*ncell];
		alpha=1.0;
		//Fin de la construccion
		
		//Importamos las posiciones y velocidades
		for(int i=0;i<npart;i++)
		{
			for(int j=0;j<3;j++)
			{
				double rr;
				MiArchivo>>rr;
				part[i].pos[j]=rr;
			}
			for(int j=0;j<3;j++)
			{
				MiArchivo>>part[i].vel[j];
			}
			boundary(i);
		}

		MiArchivo.close();
	}
}

void caja::printConfig(char namefile[30])
{
	char name[50];
	sprintf(name, "%s%s%d%s", namefile,"_",nacc,".txt");
	ofstream MiArchivo (name,ios::out);
	MiArchivo.precision(4);

	

	if (MiArchivo.is_open())
	{
		MiArchivo.setf(ios::fixed);
		MiArchivo<<"Parametros de la simulacion"<<endl;
		MiArchivo.fill('*');
		MiArchivo << setw(122) << ""<<endl;
		MiArchivo.fill(' ');
		MiArchivo<<setw(10)<<"Npart"<<setw(10)<<"phi"<<setw(10)<<"rcut"<<setw(10)<<"Dt"<<setw(10)<<"restart"<<setw(10)<<"ncorr"
			<<setw(10)<<"trelaj1"<<setw(10)<<"bx"<<setw(10)<<"by"<<setw(10)<<"bz"<<setw(10)<<"tiempo"<<setw(10)<<"paso"<<endl;
		MiArchivo.fill('-');
		MiArchivo << setw(122) << ""<<endl;
		MiArchivo.fill(' ');
		MiArchivo<<setw(10)<<npart<<setw(10)<<phi<<setw(10)<<rcut<<setw(10)<<param.Dt<<setw(10)<<param.rest<<setw(10)<<param.ncorr
			<<setw(10)<<param.trelaj1<<setw(10)<<lbox[0]<<setw(10)<<lbox[1]<<setw(10)<<lbox[2]<<setw(10)<<param.Dt*nacc<<setw(10)<<nacc<<setw(10)<<endl;
	}
	MiArchivo<<endl;
	MiArchivo<<"Posiciones y velocidades"<<endl;
	MiArchivo.fill('*');
	MiArchivo << setw(122) << ""<<endl;
	MiArchivo.fill(' ');
	MiArchivo<<setw(20)<<"rx"<<setw(20)<<"ry"<<setw(20)<<"rz"<<setw(20)<<"vx"<<setw(20)<<"vy"<<setw(20)<<"vz"<<endl;
	MiArchivo.fill('-');
	MiArchivo << setw(122) << ""<<endl;
	MiArchivo.fill(' ');
	for(int i=0;i<npart;i++)
	{
		for(int j=0;j<3;j++)
			MiArchivo<<setw(20)<<part[i].pos[j];
		for(int j=0;j<3;j++)
			MiArchivo<<setw(20)<<part[i].vel[j];
		MiArchivo<<endl;
	}

	//sprintf(aux,"%s%d%s","conf_",nacc,".txt");
	//ofstream MiArchivo (aux,ios::out);
	//if (MiArchivo.is_open())
	//{
	//	MiArchivo<<nacc<<" "<<npart<<" "<<phi<<" "<<temp<<" "<<rcut<<endl;
	//	MiArchivo<<lbox[0]<<" "<<lbox[1]<<" "<<lbox[2]<<endl;
	//	MiArchivo<<posPart;
	//	MiArchivo<<velPart;
	//	MiArchivo.close();
	//}
}

void caja::SetRandomPart(void)
{
	double pos[3];

	for(int i=0;i< npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			pos[j]=lbox[j]*double(rand())/RAND_MAX;
		}
		part[i].setPosition(pos);
		boundary(i);
	}
}

void caja::SetMeshPart(void)
{
	int np = int(ceil(pow(npart,1./3.)));
	double lx = lbox[0]/double(np);
	int n=0;
	double pos[3];
	double drift=0.5*lx;

	pos[0]=-0.5*lbox[0]+drift;
	for(int ix=0;ix<np;ix++)
	{
		pos[1]=-0.50*lbox[0]+drift;
		for(int iy=0;iy<np;iy++)
		{
			pos[2]=-0.50*lbox[0]+drift;
			for(int iz=0;iz<np;iz++)
			{
				for(int k=0;k<3;k++)
				{
					pos[k]+=(double(rand())/RAND_MAX-0.5)*0.01;
				}
				part[n].setPosition(pos);

				boundary(n);
				pos[2]+=lx;
				n++;

			}
			pos[1]+=lx;
		}
		pos[0]+=lx;
	}
	cout<<"lbox: "<<lbox[0]<<endl;
	cout<<"lx: "<<lx<<endl;
	cout<<"np: "<<np<<"  "<<int(ceil(pow(npart,1./3.)))<<endl;
}

void caja::initForces(void)
{
	for (int i=0;i<npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			part[i].force[j]=0;
		}
	}
}

void caja::followPart(int ind)
{
	cout<<"pos: ";
	for(int i=0;i<3;i++)
		cout<<part[ind].pos[i]<<"  ";
	cout<<endl;
	cout<<"vel: ";
	for(int i=0;i<3;i++)
		cout<<part[ind].vel[i]<<"  ";
	cout<<endl;
	cout<<"force: ";
	for(int i=0;i<3;i++)
		cout<<part[ind].force[i]<<"  ";
	cout<<endl;
}

void caja::linkCell(void)
{
	int mc[3], c;
	//Reset the headers, head
	for (int i=0;i<ncell3;i++)
	{
		head[i] = -1;
	}

	///Scan atoms to construct headers, head, & linked lists, lscl
	for (int i=0;i<npart;i++)
	{
		///Vector cell index to which this atom belongs
		for (int j=0;j<3;j++)
		{
			mc[j] = int((part[i].pos[j]+0.5*lbox[0])/lgrid);
			if(mc[j]>=ncell)
			{
				cout<<"mal: "<<lbox[0]<<"  "<<0.5*lbox[0]<<"  "<<part[i].pos[j]
				<<"  "<<(part[i].pos[j]+0.5*lbox[0])<<"  "<<
					(part[i].pos[j]+0.5*lbox[0])/rcut<<"  "<<ncell<<endl;
			}
		}
		///Translate the vector cell index, mc, to a scalar cell index"""
		//
		c = int(mc[0]*ncell2+mc[1]*ncell+mc[2]);
		part[i].celda=c;
		part[i].celdaxyz[0]=mc[0];
		part[i].celdaxyz[1]=mc[1];
		part[i].celdaxyz[2]=mc[2];

		///Link to the previous occupant (or EMPTY if you're the 1st)
		link[i] = head[c];
		///The last one goes to the header"""
		head[c] = i;
	}

	int ntotal=0;
	for(int i=0;i<ncell3;i++)
	{
		int cual = head[i];
		int n=0;
		while(cual!=-1)
		{
			cual=link[cual];
			n++;
		}
		ntotal+=n;
	}
	cout<<"Numero total de particulas colocadas en Mesh "<<ntotal<<endl;
}

void caja::listaVecinos(void)
{
	double r2, rc2;
	double shield = 0.5;
	rc2=(rcut+shield)*(rcut+shield);
	for (int i=0;i<npart;i++)
	{
		part[i].nvecinos=0;
	}

	for (int i=0;i<npart-1;i++)
	{
		for(int j=i+1;j<npart;j++)
		{
			r2 = distancia(i,j);
			if(r2<rc2)
			{
				part[i].lstVecinos[part[i].nvecinos]=j;
				part[i].nvecinos++;
				part[j].lstVecinos[part[j].nvecinos]=i;
				part[j].nvecinos++;


			}
		}
	}
}

void caja::initForcesKR(void)
{
	for (int i=0;i<npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			part[i].forceK[j]=0;
			part[i].forceR[j]=0;
		}
	}
}

void caja::initVelocities(double Dt)
{
	//double v0=sqrt(3.0*kB*temp);
	double v02=0.0;
	double dv[3];
	double enerkin,enerkin2;

	dv[0]=dv[1]=dv[2]=0.0;

	for (int i=0;i<npart;i++)
	{
		enerkin=0.;
		for (int j=0;j<3;j++)
		{
			part[i].vel[j]=double(double(rand())/RAND_MAX-0.5);
			enerkin+=part[i].vel[j]*part[i].vel[j];
		}
		enerkin2=0.;
		for (int j=0;j<3;j++)
		{
			part[i].vel[j]*=sqrt(3.0/enerkin);
			dv[j]+=part[i].vel[j];
			enerkin2+=part[i].vel[j]*part[i].vel[j];
		}
	}

	for(int i=0;i<3;i++)
		{dv[i]/=npart;}

	///Eliminamos el impulso inicial del sistema
	for (int i=0;i<npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			part[i].vel[j]-=dv[j];
			v02+=part[i].vel[j]*part[i].vel[j];
		}
	}

	///Transformamos las velocidades para tener la temperatura de partida deseada
	//double T = v02/(3.0*kB*npart);
	//double T = v02/(3.0*npart);
	//double factor=sqrt(1.0/T);
	//for(int i=0;i<npart;i++)
	//{
	//	for(int j=0;j<3;j++)
	//	{
	//		///Reescalamos las velocidades para tener la temperatura deseada
	//		part[i].vel[j]*=factor;
	//		///Inicializamos la posicion anterior de las particulas
	//		part[i].posOld[j]=part[i].pos[j]-part[i].vel[j]*Dt;
	//	}
	//}
	vel[0]=vel[1]=vel[2]=0.0;
	temp=temp0;
}

void caja::boundary(int ind)
{
	for(int i=0;i<3;i++)
	{
		part[ind].pos[i]-=lbox[i]*anint(part[ind].pos[i]/lbox[i]);
		//int nper = int(dabs(part[ind].pos[i])/(0.5*lbox[i])); //si esta dentro de la caja 0; si esta fuera 1
		//double anper = sign((double(nper)*lbox[i]),part[ind].pos[i]); //0 dentro de caja, 1 � -1 si esta fuera
		//Condiciones de contorno peri�dicas
		//part[ind].pos[i]-=anper;
	}
}

double caja::distancia2(int ind1,int ind2)
{
	//Calcula la distancia entre las particulas i y j
	double dis=0;
	for (int i=0; i<3;i++)
	{
		dr[i]=part[ind2].pos[i]-part[ind1].pos[i];
		drnb[i]=dr[i];
		dr[i]-=lbox[i]*anint(dr[i]/(lbox[i])); //Entero + proximo manteniendo signo 
		//dr[i] = part[ind1].pos[i]-part[ind2].pos[i];
		///Condiciones de contorno periodicas
		//dr[i]-=int(dr[i]/(0.5*lbox[i]))*lbox[i];
		dis+=dr[i]*dr[i];
		//cout<<dr[i]<<" ";
	}
	//cout<<dis<<endl;
	return	 dis;
}

double caja::distancia(int ind1,int ind2)
{
	return	 sqrt(distancia2( ind1,ind2));
}

void caja::printXYZfile(char namefile[30])
{
	char name[50];
	sprintf(name, "%s%s", namefile,".xyz");
	ofstream MiArchivo (name,ios::out|ios::app);
	MiArchivo.precision(5);
	if (MiArchivo.is_open())
	{
		MiArchivo << npart<<endl;
		MiArchivo << lbox[0]<<"  "<<nacc<<endl;
		for(int i=0;i<npart;i++)
		{
			MiArchivo<<setw(5)<<"O "<<setw(8)<<part[i].pos[0]<<"  "<<setw(8)<<part[i].pos[1]<<"  "<<setw(8)<<part[i].pos[2]<<endl;
		}
		MiArchivo.close();
	}
}

void caja::printShearStressTensor(char namefile[30])
{
	char name[50];
	sprintf(name, "%s", namefile);
	ofstream MiArchivo (name,ios::out|ios::app);
	if (MiArchivo.is_open())
	{
		MiArchivo << nacc*param.Dt<<"  "<<(svxz+spxz+svyz+spyz)/2.0<<"  "<<svzz+spzz<<endl;
		MiArchivo.close();
	}
}

void caja::iniVelocidades2(void)
{
	double v2[3];
	double enerkin, f;
	ekin=0.;
	for(int i=0;i<npart;i++)
	{
		enerkin=0.;
		for(int j=0;j<3;j++)
		{
			part[i].vel[j] =3.0*(double(rand())/RAND_MAX-0.5);
			v2[j]=(part[i].vel[j]*part[i].vel[j]);
			enerkin+=v2[j];
		}
		f=3.0/enerkin;
		for(int j=0;j<3;j++)
		{
			part[i].vel[j]*=sqrt(f);
			v2[j]=(part[i].vel[j]*part[i].vel[j]);
			ekin+=v2[j];
		}
	}

	vel[0]=vel[1]=vel[2]=0.0;
	ekin*=0.5/double(npart);
	temp=2.0/3.0*ekin;
//	ekin=1.5;
}

void caja::SetMeshPart2(void)
{
	int nmax3=anint(pow(npart,1./3.));
	int ipart=0;
	double lg=lbox[0]/(pow(double(npart),(1./3.)));
	double alea;

	for(int i=1;i<=nmax3;i++)
	{
		for(int j=1;j<=nmax3;j++)
		{
			for(int k=1;k<=nmax3;k++)
			{
				alea=(double(rand())/RAND_MAX-0.5);
				part[ipart].pos[0]=-lbox[0]/2.+((2.0*double(i)-1.0)/2.)*lg+lg*(1./4.)*alea;
				alea=(double(rand())/RAND_MAX-0.5);
				part[ipart].pos[1]=-lbox[1]/2.+((2.0*double(j)-1.0)/2.)*lg+lg*(1./4.)*alea;
				alea=(double(rand())/RAND_MAX-0.5);
				part[ipart].pos[2]=-lbox[2]/2.+((2.0*double(k)-1.0)/2.)*lg+lg*(1./4.)*alea;
				ipart++;
			}
		}
	}
}